The X-Philes (2nd Revision)/The X-Philes Number 1 (1995).iso/xphiles/hp48_1/stoer/comp.sys.handhelds_2970_000000.msg

Internet Message Format | 1995-03-31 | 7KB

From: [e-mail] (Matt Austern)
Subject: Bulirsch-Stoer integration on the HP 48  (Part 1 of 2)
Summary: Documentation, description, and examples.
Keywords: Bulirsch-Stoer, differential equations
Date: 18 Dec 90 18:13:28 GMT
Organization: Lawrence Berkeley Laboratory (theoretical physics group)

I have written a program for the hp 48 which solves systems of
differential equations using the Bulirsch-Stoer algorithm.  This
posting is a description of the program, and the next posting is the code
itself.

SUMMARY OF THE METHOD

If you look at the code, it will immediately be apparent that this
algorithm is immensely more complicated than, say, Runge-Kutta.
Although this complexity makes a single step rather slow, the
advantage of Bulirsch-Stoer is that it is often possible to choose
extremely large step sizes while maintaining reasonable accuracy.

The Bulirsch-Stoer algorithm is described in Numerical Recipes, by
Press, Flannery, Teukolsky, and Vettering.  (If you do any numerical
work and you don't already have a copy of that book, by the way, you
should get one.)  I'll give a brief summary, just so that this
program won't be completely a black box.

If you have a system of equations of the form y' = f(x, y), where y
and y' are in general vectors, and the initial conditions (x0, y0),
then the problem is to find y(x1) for some x1.  (Normally, we require
x1 - x0 to be small; here, we make no such requirement.)  Define H =
x1 - x0.

We divide H into n subintervals, and then use n iterations of a
simple first-order method (specifically, the modified midpoint
method) to find y(x1).  Define this value to be y(x1, h), where
h = H/n.

The Bulirsch-Stoer algorithm is to compute y(x1, h) for several
different values of h, and then extrapolate down to h=0.  The
extrapolation algorithm provides an estimate of its accuracy; when
this accuracy is good enough, we return the answer.  Convergence is
somewhat quicker than might be expected; it turns out that
y(x1, h) contains only even powers of h, so the extrapolation is
actually in h^2.

This algorithm can occasionally fail.  The rational function
extrapolation might encounter a pole, in which case there's nothing
to do but quit.  The program checks for this failure, so that at
least there won't be a mysterious blowup.  (This is rare, but it
happened to me once.  The solution is to use a different
extrapolation procedure in those cases.  If I get energetic, I may
implement polynomial extrapolation for the 48sx as a fallback.)  It
is also conceivable that the step size will have to be reduced so far
that x+h will be indistinguishable from x.  I have never seen this
happen (my guess is that it would only happen when you're close to
some singular point in the solution of the differential equation, or
when you have a system of equations which vary on very different
scales), but there is code to check for that too.

THE PROGRAM:

There are actually four programs here: EXTRAP, NDESTEP, DESTEP1, and
DE.  (There are also two objects that are intended for internal use.
They are named in lowercase.)

EXTRAP is the extrapolation routine.  It takes two arguments, both
lists.  Level 2 contains a list of x values, and level 1 a list of y
values.  EXTRAP returns two values: in level 2 it returns y(0), and
in level 1 an estimate of the error in that prediction.  (EXTRAP,
incidentally, is where most of the time gets spent.  I've tried to
make the code in its inner loop efficient, but I'd welcome any
suggestions for speeding it up.)

DESTEP1 solves a single first-order differential equation.  It takes
five arguments.  In order, level 5 to level 1, they are: ydot,
tolerance, stepsize, x, y.  ydot is a function: it takes x and y (on
the stack, in that order) and returns y'(x, y).  Tolerance is the
required accuracy (i.e., if Delta-y is the error, then we demand
Delta-y / y < tolerance), and stepsize, x, and y are
self-explanatory.  DESTEP1 returns the same information, to make it
convenient to take another step.  ydot and tolerance are left
unchanged on the stack; the new step size is that recommended by the
algorithm (you don't have to follow the recommendation, of course!);
and x and y are the new values.

(Note: if you are too greedy with your step size and your tolerance,
DESTEP1 might have to choose a smaller step than you asked for.  If
so, x will be incremented by the step actually used.  This is
slightly unusual, though: Bulirsch-Stoer can handle most reasonable
requests.)

NDESTEP is just the same, except that y is a vector, and ydot must
take a vectorial y as an argument and return a vectorial ydot as a
result.  (Most of the code in NDESTEP, actually, is identical to that
in DESTEP1.  Merging them would be quite easy if you want to conserve
memory.)

DE is intended for interactive use; it's essentially just a
cosmetic shell for NDESTEP and DESTEP1.  It takes and returns only
four arguments: tolerance is taken from the display format.  (In STD
format, it uses a tolerance of 0.0001.)  It chooses to call NDESTEP
or DESTEP1, depending on the type of argument given.  It returns x
and y as tagged objects; if the user provided tags then those tags
are preserved, otherwise it uses the rather unimaginative labels "x"
and "y".  Finally, NDESTEP and DESTEP1 use user flags; DE saves and
restores the old values.


EXAMPLE PROBLEM:

Solve the equation x^2 y''(x)  +  x y'(x)  +  x^2 y(x)  =  0, subject
to the initial conditions y(0) = 1, y'(0) = 0.  Specifically, find y(5).

Solution:

This is equivalent to the system
    y1'  =  - y1/x  -  y
    y'   =  y1,
with y(0) = 1, y'(0) = 0.

Put the following four entries on the stack:
\<< OBJ\-> DROP \-> x y1 y
    \<< IF x 0 == THEN y NEG ELSE y1 x / y + NEG END
        y1 2 \->ARRY
    \>>
\>>
5
0
[0 1]

Then, with the display mode set to 3 FIX, hit DE.  After 67 seconds,
I get the result that y(5) = -0.178, and y'(5) = 0.328.

In fact, the solution to this equation is a Bessel function, J0(x).
The tabulated answer is J0(5) = -0.177597, and J0'(5) = -J1(5) = 0.327579.

67 seconds is a long time, but we were able to span the entire
interval from 0 to 5 in a single step, and still get three digits of
accuracy.


Enjoy!  There are obviously better tools than the HP 48 for serious
numerical calculations, but it's very convenient to have a quick way
of playing with a differential equation.
--
Matthew Austern    [e-mail]     Proverbs for paranoids, 3: If
[phone]     [e-mail]    they can get you asking the wrong
                   [e-mail]        questions, they don't have to worry
                                          about answers.